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The Joint Use of the Ordinary and Extraordinary Virtual 
Height Curves in Determining Ionospheric Layer Profiles 
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An extension of Budden's matrix method for determining ionospheric layer profiles is 
described. When analyzing vertical incidence ionograms by the matrix method, it is usual 
to interpret the virtual height curve for the ordinary mode only. Errors then arise from 
the presence in the lower ionosphere of low-density ionization for which the plasma frequency 
is less than the lowest frequency observed. In the proposed extension of the method, such 
errors are reduced by making use of the extraordinary virtual height curve as well as of the 
ordinary. 



1. Introduction 

In 1954 Budden [1] 2 introduced the matrix 
method for analyzing the records of vertical incidence 
ionospheric soundings, and since then it has come 
into wide use. The records, which are known as 
"ionograms," contain curves of the virtual reflection 
heights of the ordinary and extraordinary modes 
as functions of frequency ('Virtual height curves"), 
and the object of the analysis is to determine, from 
one or the other of these curves, how the ionization 
in the underside of the ionosphere is distributed in 
height. The difficulty lies in the integral equation 
that relates this distribution to the virtual height 
curve; in general, it is too complex to be solved 
analytically. The principle of Budden's method is 
to approximate the integral by a sum of discrete 
values of the integrand, thus converting the integral 
equation into a matrix equation, which may be 
solved by inverting the matrix. 

In its customary form, the method makes use 
only of the ordinary virtual height curve. The 
ordinates of this curve are measured from the 
ionogram at selected frequencies, and then this set 
of ordinates is multiplied by a matrix to produce a 
set of the corresponding ordinates of the layer 
profile. The matrix multiplication is a simple 
repetitive process and so lends itself readily to 
performance by a digital computer. As a result, 
the matrix method is now widely favored for the 
routine analysis of ionospheric data. 

Unfortunately, however, two rather restrictive 
conditions have to be met if the method is to work 
properly: First, the original virtual height curve 
should extend down to zero frequency; second, the 
ionosphere should consist of a single layer only, so 
that a given density of ionization occurs at just 
one height on its underside. Only under these 
conditions does the integral equation define a unique 
distribution. 
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Usually neither condition is met in practice. 
Most ionospheric sounders do not operate below 
about 1 Mc, while the ionosphere often appears 
to be stratified into several distinct layers. Then 
there are many possible distributions that could 
have given rise to the observed portion of the virtual 
height curve. In fact, these uncertainties are the 
main sources of error in the analysis of ionograms 
by the matrix method. 

This paper deals only with the problem of how 
to reduce the errors that arise from ignorance of 
the virtual heights at low frequencies. Here the 
ionosphere still is assumed to consist of a single 
layer. The analysis appears capable of being 
extended to cover the case of multiple layers with 
intervening "valleys," but this development is not 
considered. 

When the virtual height curve is missing at the 
low-frequency end, the usual practice is to assume 
that the unknown virtual heights are constant and 
equal to the value at the lowest frequency observed; 
this assumption implies that there is no ionization 
below the level of reflection for that frequency. 
If such "low-density ionization" is present in fact, 
then the group delay that it produces in the observed 
frequency range is attributed wrongly to time spent 
by the waves in traveling up to the ionosphere, so 
that the true heights are overestimated. 

So long as attention is confined to a single virtual 
height curve, there is no way of detecting the low- 
density ionization or of isolating its effects in the 
observed frequency range. But these effects are 
different for the two modes of propagation; the 
group delays differ in amount and vary with fre- 
quency in different ways. This fact suggests that 
the low-density ionization could be studied by 
measuring both the ordinary and extraordinary 
virtual height curves and combining the data that 
they contain. In the present paper, extensions of 
the Budden matrix method arc developed to do this. 

The contents of the paper are arranged as follows: 
In section 2 terms are defined, and the present method 
for the interpretation of a single virtual height curve 
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is reviewed briefly; the work reviewed is chiefly 
that of Budden [1], Rydbeck [2], and Schmerling [3], 
Section 3, which contains the substance of the paper, 
shows how the matrix method may be extended 
to make joint use of both virtual height curves; the 
errors that may arise in the inferred distribution 
of ionization, due to random errors in measuring 
the virtual heights, are examined also. The analysis 
is illustrated in section 4 by a worked example, in 
which an artificial ionogram, computed for a known 
distribution of ionization, is interpreted by the 
new procedures. Finally, the whole analysis is 
discussed and criticized in section 5, and the con- 
clusions of the paper are summarized in section 6. 
Certain previous authors also have considered 
using the two virtual height curves in conjunction. 
The suggestion was made originally by Jackson [4], 
but his approach, which is one of trial and error, 
is not suitable for extensive use. Thomas [5], 
however, in a recent survey paper, has outlined some 
work by Titheridge [6], who appears to have adapted 
Jackson's method into a more convenient form. 
Titheridge ? s work, which has not been published at 
the time of writing this paper, ma} r anticipate some 
of the results presented here. 

2. Interpretation of a Single Virtual Height 
Curve 

2.1. The Integral Equations 

In this section the integral equations are derived 
that relate the virtual height curves for the ordinary 
and extraordinary modes to the distribution of 
ionization. 

At vertical incidence, the virtual height of reflec- 
tion hf is given as a function of the wave frequency 
/ by the integral 



h'(j)=j\'(f,N)dz, 



(1) 



where z is height in general, z R is the true height of 
reflection for the given frequency, and \i! is the 
group refractive index, which depends both on the 
frequency and on the local electron density N. 

In the limit of very low wave frequencies, the 
virtual and true heights of reflection for both modes 
all tend to a single limiting value z , the height of the 
base of the ionosphere. From ground level up to 
this height the waves travel at the speed of light, 
which is to say that their group refractive indices are 
unity. 3 Thus (1) can be written as 



h'(f)=2o+ ( ZR n'(f,N)dz 



(2) 



3 Actually /*' is slightly greater than unity even below the ionosphere due to 
the refraction of the air, but the departure from unity is negligible in the present 
context. 



where the integral now embraces only the range of 
height in which free electrons are present. 

The general expression for the virtual height is com- 
pleted by introducing the relation between the density 
of electrons and the height. However, instead of the 
density itself, it is more convenient to work with the 
"plasma frequency" J N . which is proportional to its 
square root. The function z(J N ), in which plasma 
frequency is viewed as the independent and height as 
the dependent variable, will be referred to here as the 
"layer profile." The profile is a multivalued function, 
even for a single layer, since a given value of j N , if it 
occurs at all, must do so at least twice, once on the 
upper surface of the layer and once on its underside. 
The relationship between the true heights of reflection 
and the layer profile may be understood by consider- 
ing the condition for reflection. The condition for an 
ionized medium to reflect radio waves is that its 
plasma frequency should exceed a certain critical 
value that depends on the frequency of the waves 
and also on the mode of propagation. This value will 
be written as /«(/). It follows that waves of a certain 
frequency and mode, propagated up into the iono- 
sphere from below, are reflected at the first level that 
they encounter at which f N is equal to the appropriate 
value of f R . Obviously, this is the lowest such level. 
Hence, the true height of reflection is given by 



z R = minimum value of z{j R ). 



(3) 



To summarize, the true height of reflection is a 
particular value of z, determined by f R , which in 
turn is a function of /; this function is different for 
the two modes of propagation. 

Since all reflection takes place from the underside 
of the ionosphere, only the profile of this part needs 
to be considered. Further, since the ionosphere has 
been assumed to consist of a single layer, the profile 
of its underside is a mono tonic function. In this 
case, f N can be taken as the variable of integration. 
Also, the group refractive index // will be expressed 
as a function of f N instead of N. With these changes, 
the integral equation becomes 



h' (/) = Zo+ JV {J, f«) (J£) df N . 



(4) 



Here the equation is still quite general, in that it 
applies equally to both modes. The particular 
equations for the ordinary and extraordinary modes 
are to be obtained by considering their respective 
reflection conditions, that is to say, the different 
relationships between f R and /. These relationships 
are graphed in figure 1. 

The ordinary mode behaves very simply; the critical 
plasma frequency for reflection is equal to the wave 
frequency. Thus, 



^R=f (ordinary mode) 



(5) 
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Figure 1. The conditions for reflection of the ordinary (0) 
and extraordinary (X) magneto-ionic modes at vertical 
incidence. 

Abscissa: wave frequency (/). Ordinate: plasma frequency at the level of 
reflection (fa). Both frequencies are normalized with respect to the gyrofre- 
quency (fit). 



so the integral equation for the ordinary mode is 

*' (/) =So+ J* V (/, f K ) (j£) dU (6) 

From here on, the quantities h, h' ', and //, when 
written without subscripts, will refer to the ordinary 
mode only. The corresponding quantities for the 
extraordinary mode will bear the subscript "x." 

The extraordinary reflection condition is more com- 
plex. The relationship between f R and /is now 

h = { /(/± Jh) }^ (extraordinary mode) (7) 



=/, 



say 



(8) 



where in is the gyrofrequency. The positive sign 
applies if /</# (L.H. curve in fig. 1), while the 
negative sign applies in the more usual situation, 
where />/# (R.H. curve). 

Note that, as the wave frequency approaches the 
gyrofrequency from below, f R tends to the value -yj2 
f H . On the other hand, as the gyrofrequency 7 is 
approached from above, j R becomes vanishingly small, 
so that the extraordinary reflection level should de- 
scend to the base of the ionosphere; actually, such 
echoes are never detected, because they are heavily 
absorbed . 

The integral equation for the extraordinary mode 
can be written now as 



K (/) = 2 +f/ x (JJn) (Jr)#* 



(9) 



This equation differs slightly from the correspond- 
ing equation for the ordinary mode, in that the upper 

limit of the integral (/) is not the same as the inde- 
pendent variable (f) on the left-hand side; also, to 
each value of/ below the limit -\f2f H there correspond 
two possible values of/ (see fig. 1), one of which is 
less than the gyrofrequency, while the other is 
greater. However, only the upper range of / will 
be considered here, because most records of iono- 
spheric virtual height do not extend below the gyro- 
frequency. With this restriction, there is a straight- 
forward correspondence between / and /, and the 
virtual height may be regarded equally well as a 
function of either variable. When it is regarded as 
a function of/, the extraordinary virtual height will 
be represented by the symbol h; thus 



h(f)=K(f). 



(10) 



where the symbol / has been adopted to denote the 
value of f R for an extraordinary wave of frequency/. 



To find the virtual height for a specified value of 
/, it is necessary to know the corresponding value of 
/. This frequency is given by the positive solution 
of the quadratic eq (7), and it is 

f=ifH{{i+Mf/f H ) 2 } i Ti\, (li) 

where the positive sign now yields the range of 
frequencies above f H , which is the range to be 
considered. 

In these terms, the integral equation for the extraor- 
dinary mode is 

kf)=Zo+j^' x (f,U) (|r)#*- ( l2 ) 

This equation is similar to (6), the only remaining 
difference being that, in calculating y! x , it is neces- 
sary first to obtain from (11) the value of/ that 
corresponds to the given value of/. 

The integral eqs (6) and (12) relate the virtual 
height curves to the slope of the layer profile. Were 
it possible to invert either equation so as to obtain 
a solution for the slope, then the profile itself could 
be recovered by a straightforward integration. 

2.2. The Matrix Equations 

The integral equations cannot be solved analyt- 
ically in the general case, because the group refrac- 
tive indices are complicated functions of / and j N . 
Therefore, they are solved by converting them into 
matrix equations as will be described now. 

Consider the eq (6) for the ordinary mode. The 
first step is to replace the continuous functions h r 
and z by their values at a number of discrete "sample 
frequencies. " These frequencies will be written, in 
order of ascending magnitude, as/ , /i,/ 2 . . . etc., 
or generally as j m where the index m is zero or a 
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positive integer; J is zero frequency, and the rest 
are assumed to be spaced closely, but not necessarily 
evenly, through the full frequency range. The fol- 
lowing notation will be used also : 



&Jm — Jm Jm — 1 


(m^l) 


(13) 


K=h'(f n ) 




(14) 


Zm — 2\Jm) 




(15) 




(ro^l) 
(m=0). 


(16) 



The quantities Az m are increments of height on the 
layer profile; note the definition for m=0. 

The next step in the approximation is to assume- 
that the slope of the layer profile is effectively 
constant over each elementary frequency interval 
Aim, i-e., 



dz ^Az m 
dj N ~Aj m 



{jm-aJN^fm). 



Then the complete integral can be split up into a 
series of discrete sums, as follows: 

m = \L\J m JJm-\ 



^Zl M nm Az m , 

m = Q 



say 



(17) 



where Az has been written in place of z , and 
r \ (m=0) 



M nm =« 



4/, 




T \ /'" / (JuJn) dj N (lZmZn) (18) 

JmJ Jm—1 



Now suppose that in (17), n is set equal in turn to 
0, 1, 2, etc., up to the value that corresponds to the 
highest sample frequency. The result is an array of 
simultaneous equations, which may be solved for 
the unknowns Az m by the simple algebraic process 
that is described below in section 2.4. Originally, 
however, Budden [1] adopted an alternative ap- 
proach, which is to regard the complete array as a 
single "matrix equation ;" this approach has pro- 
vided a more convenient basis for the subsequent 
extensions of Budden's methods. In matrix termi- 
nology, the set of quantities ti n and also the set of 
Az m , both constitute "vectors/' while the coefficients 
M nm form a "lower triangular matrix." This differs 
from the corresponding matrix used by Budden 
(1955) and Schmerling (1957) by including an extra 
column with m=0, in which all the coefficients are 
unity, and an extra row with n=0, in which all the 
coefficients except the first are zero. Each of the 
remaining nonzero coefficients is the average of the 



values of group refractive mdex assumed by a wave 
of frequency f n in traveling from the level where 
j=f m -\ up to that where f N =f m . 

If n is understood to run through its full range of 
values, then eq (17) states that the vector of the 
h' n is obtained by multiplying the vector of the Az m 
by the matrix of the coefficients M. nm . Thus (17) is 
the required matrix equation for the ordinary mode. 

The integral equation for the extraordinary mode 
can be treated likewise, and the result is the matrix 
equation 



An^S M nm Az n 

71 = 



(19) 



where 



f\ 



^-nm — * 



(m=0) 



1 I f m i / _f 

—ST ~ HxUn 

AfJ f ^-i 



, In) df N 







(l^m^n) (20) 
(m>7i) 



and the other new quantities are defined like their 
ordinary counterparts (cf. eqs (13), (14)). 

From here on, the complete matrices will be rep- 
resented by symbols in square brackets, while the 
same symbols without brackets will refer to the 
individual elements. Thus the group refractive index 

matrices will be written as [M nm ] and [M ww ]. Their 

elements M nm and M nm may be computed from the 
formulas given in appendix A. 

2.3. Solution by Matrix Inversion 

The matrix eqs (17) and (19) for the ordinary and 
extraordinary modes, respectively, may be solved 
by inverting the matrices. Thus, for the ordinary 
mode, the solution is 



Az n 



w = 



(21) 



where the matrix [N nm ] is the inverse of [M nm ] . Sim- 
ilarly, for the extraordinary mode, 



n „ „ 



I\Z n / , IV nm fom 

m = Q 



(22) 



where [N nm ] is the inverse of [M nm ], Because the 

matrices [M nm ] and [M nm ] are lower triangular, so 
also are their inverses. 

The increments of height may be summed now to 
give the heights themselves: 



Zn=12 AZ 
m=0 



(23) 



This summation is analogous to the final integration 
that would be performed if the original integral equa- 
tions, that involve the slope dz/dj N , could be solved 
analytically (see sec. 2.1). 



114 



2.4. Solution Row- by-Row 

As an alternative to inverting the matrices, the 
increments of height on the layer profile may be 
found by solving the matrix equation row-by-row 
[7]. For instance, consider the eq (17) for the 
ordinary mode. It can be renmuiged to read 

Az "— XT"! *^""S M nn Az n \ . (24) 

There is one such equation for each value of n. The 
complete set may be solved by starting from n=0 
and working upwards. Thus, 



1 

~M o 



A^o^nr-Uo}' 



Azi^jf-IK— M 10 Az } 



AZ2 ~W { A 2-M 21 A3 1 -M 20 A^ }^ 

and so on. In effect, the method consists of isolat- 
ing, at each stage, the group delay produced in the 
unknown uppermost increment of height, by sub- 
tracting from the observed virtual height the group 
delay produced in the known height increments at 
lower levels. 

The eq (19) for the extraordinary mode can be 
solved similarly. Matrix equations can be solved in 
this way only when, as in the present case, the 
matrices are triangular. 



3. Joint Use of Both Virtual Height Curves 

3.1. General Remarks 

In the previous section, each of the matrix eqs 
(17) and (19) represented a set of linear simultaneous 
equations in which the num.ber of data (ordinates of 
the virtual height curve) was equal to the number 
of unknowns (increments of height on the layer pro- 
file). Hence each set had an exact solution, so that 
the complete virtual height curve for either mode 
was sufficient by itself to determine the profile. 
But when some of the curve is missing, as when the 
records are unavailable at low frequencies, then the 
remaining high frequency data are insufficient to 
determine the profile without making special assump- 
tions (see sec. 1). 

However, the situation is not hopeless, since the 
data for the two magneto-ionic modes are mutually 
independent. Suppose, therefore, that in sampling 
the virtual height curves to obtain these data, the 
sample frequencies are chosen so that the unknowns 
are the same in the two sets of equations. Then, if 



both sets are combined, the number of data is in- 
creased while the number of unknowns remains con- 
stant. Hence, the total data may now be enough 
for the complete set of equations to be soluble. 

To make the unknowns the same, the ordinary 
and extraordinary virtual height curves must be 
sampled at frequencies which correspond to the same 
true heights of reflection. Considering one such 
pair of frequencies, this is to say that the values of 
j R must be the same for the two modes. From the 
expressions for j R ((5) and (7)), this requirement is 
seen to be met if the sample frequency j n for the or- 
dinary mode is made equal to j n for the extraordi- 
nary; the corresponding sample frequency /„ for the 
extraordinary may be found from (11). Thus the 
corresponding values of j n are not equal for the two 
modes; moreover, if the sample frequencies are 
spaced evenly for the ordinary mode, then those for 
the extraordinary mode will be spaced unevenly. 
From here on, the two sets of sample frequencies will 
be assumed to be related in this way. 

The question now arises, under what condition 
can the layer profile be determined fully by joint 
use of both virtual height curves? Suppose, to 
simplify the discussion, that all the data come in 
pairs; sample values of ordinary virtual height for 
which no corresponding extraordinary virtual heights 
are measurable, or vice versa, will be ignored. 
Suppose also that the range of index n for the two 
sets of data is r^ri^r+s-l. Then the data com- 
prise s values of h' n , together with s values of h n , 
a total of 2s separate items. The number of the 
unknowns Az m is r-\-s. Evidently the condition 4 
under which the data determine the unknowns fully 
is that s^tr. 

If the positions of the pairs of sample frequencies 
can be chosen arbitrarily, it is easy to arrange for 
a surplus of data. If, on the other hand, it is decided 
to space the frequencies evenly for one or other of 
the modes, say the ordinary, then the condition 
entails that the frequency range over which the 
ordinary virtual height is measurable should be 
wider than the unobserved range of low frequencies. 
Since, for most ionosondes, the low frequency limit 
lies below 2 Mc, while the ordinary penetration 
frequency for the F2 layer usually exceeds 4 Mc, 
this requirement is likely to be met most of the time. 

Assume, then, that there are enough data to 
determine the unknowns. Two distinct cases are 
to be considered: In the first, s is equal to r, so that 
the data determine the unknowns exactly; in the 
second, 5 is greater than r, and the unknowns are 
over-determined. The exactly-determined case is 
treated in section 3.2, and the overdetermined case 
in section 3.3. 



4 Strictly speaking, this condition is necessary but not sufficient, for, even when 
it is satisfied, the determinant of the set of equations still might vanish. In 
practice, the determinant is unlikely to vanish exactly, though it might become 
very small. Then quite small errors in the data would produce large errors in 
the inferred values of the unknowns, a state of affairs that is to be avoided. For 
a discussion of the effects of random errors, see section 3.5. 
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3.2. The Exactly-Determined Case 

In the exactly-determined case, the index n for 

the data {h' n and h n ) runs from r to 2r-l . The matrix 
eqs (17) and (19) together represent a set of 2r 
simultaneous equations, which are to be solved for 
the 2r unknown values of the Az m . 

It is convenient to split the solution into two 
steps: First, the equations are solved for the incre- 
ments of height in the unobserved frequency range 
(0^m<r); after which, the remaining increments 
in the observed range (r^m<2r-l) are found by 
existing methods. 

In finding the Az m for the range 0^m<>, the 
values outside this range first must be eliminated. 
This is done most simply by starting from the in- 
verse eqs (21) and (22). Because the sample 
frequencies for the ordinary and extraordinary 
virtual height curves correspond to the same true 
heights of reflection, these two expressions must 
give the same values for the Az n to the limits of 
accuracy of the matrix formulas. Therefore, they 
can be equated to yield the relation that holds 
between the two virtual height curves as a conse- 
quence of their being derived from the same layer 
profile : 



2—1 X-^nrnf^n 
m=0 



-N h' ] ^0 



(25) 



When n^r, (25) provides a relation between the 
observed and the unobserved virtual heights. This 
fact can be made plain by splitting the sum into 
separate parts for the two frequency ranges: 

n „ „ r— 1 „ ~ 

HZiNnJl'm—NnmK} ~ 2 {N nm h m — N nm fQ (n > r) . 
m = r m=0 

(26) 

The next step is to substitute for the unobserved 
virtual heights and reintroduce the corresponding 
Az m by use of (17) and (19): 

X! {N nm h' m -N nm hm} ^S i N n J SMni Az, ) 

m=r m=0 V. \l=0 / 



N nm (jtM m , AzM- (n>r) 



r— 1 m „ „ 

■1^1^{N nm M ml -N nm M ml }Az l (n>r) (27) 

m = 1=0 



=* X) S {N nm M m i-N nm M ml } As t (n> r) (28) 

1=0 m = l 

The reversal of the order of summation on the 
right-hand side, between (27) and (28), may be 
justified by pointing out that, in the form of (27) 
a given m occurs in combination with all the Z's 
from to m; hence, a given I must occur in combina- 
tion with all the ra's from / to r-1, as indicated in 



(28). Or again, it may be checked simply by 
writing down the triangular array of the possible 
combinations of I and m. 

Now (28) may be written in the concise form 



where 



S n ^J^P nl Az l (n>r) (29) 

1=0 



Sn=T>{N nm K' m -R nm l m } (n>r) (30) 



and 



Pm=lZ{N nm M nl -N nm M ml } (n>r). (31) 

m=l 



The matrix eq (29) relates the A^z for the unob- 
served frequency range to the quantities S n , which 
are formed from pairs of corresponding values of 
ordinary and extraordinary virtual height in the 
observed range. Its solution will be written as 



AzfiJ^QiJSn (l<r-l), 

n=r 



(32) 



where the matrices [Qi n ] and [P n i] are mutually in- 
verse. 

Having thus arrived at the height increments Az x 
for /<r-l, the remaining values can be found by 
solving one of the original matrix equations row-by- 
row (see sec. 2.4). Use of the virtual height data for 
either mode should lead to identical results, even 
when the data contain errors of measurement, be- 
cause the equations for the two modes together form 
an exactly-determined set. 

3.3. The Overdetermined Case 

In the overdetermined case, where s>r, any subset 
of r+s selections from the 2s items of data would 
suffice, in principle, to determine the unknowns com- 
pletely. If the data were free from error, it would not 
matter which subset was chosen. But when the data 
contain errors of measurement, some subsets will 
give better estimates than others. A still better esti- 
mate might be obtained by taking an appropriately 
weighted mean of the estimates from all the possible 
subsets, in which case each unknown would be given 
as a linear combination of all the data. The choice of 
the best combination would be governed by the 
statistics of the errors. Here it will be assumed, for 
convenience, that the errors of measurement of the 
virtual heights are all random, mutually uncorre- 
lated, and distributed normally with equal standard 
deviations. Then the best estimate of the Az x is that 
given by the "Method of Least Squares" (see, for 
instance, Whittaker and Robinson [8]), and it is the 
solution of the matrix equation 



B m ^ X) U m lAZ h 



(33) 
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where 



r+s-l „ „ 

R m = 2 {M nm h' n +M nm h„} 



(34) 



and [U m i] is a square symmetric matrix, the elements 
of which are given by the expression 



U m 



r+s-l « „ 

= S {M nm M nl +M nm M n 



(35) 



These equations are derived in appendix B. 

The matrix eq (33) represents a new set of r+s 
linear simultaneous equations, equal in number to 
the unknowns; each of the new equations is a linear 
combination of the 2s original equations, represented 
by (17) and (19). The solution of the complete set 
will be written, in matrix form, as 



r+s-l 



Az^ S V lm R n 



m=0 



(36) 



where the matrix [V tm ] is the inverse of [Ui m ] and like- 
wise is square symmetric. 

In eqs (33) to (36) , the index r would be a constant, 
fixed by the properties of the recording equipment, 
but the index s would vary from one set of data to 
another because of the variations in the penetration 
frequencies of the F2 layer. Equation (36) is the 
solution of (33) just so long as s>r. The equations 
do not contain this restriction explicitly, but when 
s<> the determinant of [U m i] vanishes so that the 
matrix cannot be inverted. Equation (36) is indeed 
applicable to the exactly-determined case when s=r, 
but for this special case the equations developed in 
section 3.2 give the same answer more simply; at 
least, the matrix that has to be inverted is smaller. 

3.4. Some Practical Considerations 

The variability of the index s makes the full least- 
squares solution inconvenient for practical use be- 
cause each value of s requires a different matrix 
Vim to be computed. In practice, it is best to fix 
s at some particular value, so that Vi m need be com- 
puted only once. The chosen value of s should be 
sufficiently high for the unknowns to be thoroughly 
overdetermined, yet sufficiently low for the upper- 
most sample frequency (/ r+s _ x ) never to exceed 
f F 2 . The use of this particular matrix [V im ] yields 
the values of the Az x for ^ I ^ r-\-s— 1 . When data 
are available at frequencies above jfr+t-i, the remain- 
ing Azi may be obtained by continuing the row-by- 
row solution upwards. If the data contain random 
errors, however, the results would be different 
according as the ordinary or extraordinary virtual 
heights were used in the solution. Better results 
would be obtained by combining both sets of virtual 
height data following the Method of Least Squares. 
The least-squares value for each of these Azi(l>r-{-s) 
is found by applying the procedure of section 2.4 



to the equation 



MR 



■S M lm Az m V 



m=0 



M, 



^r- < A,-S M lm Az m V (37) 



which replaces (24). 

Probably not much accuracy is lost by fixing the 
value of s in the way suggested, which entails that 
only data at frequencies between f r and f r + s -i are 
used to determine the distribution of the low- 
density ionization. For this ionization has its 
greatest effects on the observed virtual heights at 
frequencies just above f r , and hence these frequencies 
must provide the data that are most useful in de- 
termining its distribution. 

Furthermore, the virtual height data at the higher 
frequencies may be systematically unreliable as a 
result of "lateral deviation. " This effect is discussed 
more fully in section 5. Here it is mentioned only 
as a further justification for omitting the data at 
frequencies above f T + s -i initially, when using the 
matrix method for the overdetermined case. 



3.5. The Effect of Random Errors 

If the measurements of the virtual heights con- 
tain random errors, these errors produce fluctuations 
in the inferred values of the Az h The standard 
deviation of these fluctuations will be calculated now 
for the total range of I. The calculation will be made 
for the overdetermined case assuming that the 
Azi have been obtained by the full least-squares 
procedure ((33) to (36)). As before, the errors in the 
measurements will all be assumed to have the same 
standard deviation, and this will be called o\ 

Now the change that occurs in Az h as a result of 
a small error 5h' n in the measurement of the ordi- 
nary virtual height h' n , is found by differentiating 
(36) and (34): 



r+s-l /?\J? \ 



TO = 

r+s-l 
ra=0 



(38) 



An error in the measurement of the extraordinary 

virtual height h n has a similar effect. The total 
standard deviation of Az h due to all the errors, is 
found by squaring and adding these two contribu- 
tions, summing over all n, averaging, and then 
taking the square root. The results can be written 
as: 



Standard deviation of Azi—<rp h 



(39) 



117 



where 



Pi- 



r+s-l frr+s-l "12 rr+s-1 „ -12"N 



(40) 



The "error coefficient 7 ' p h which is the standard 
deviation of Az x per unit standard deviation of the 
original measurements, is a measure of the suscepti- 
bility of the estimate of Az t to disturbance by ran- 
dom errors of measurement; the smaller is p h the 
more reliable is the estimate. 

For the exactly -determined case, the error coeffi- 
cients may be obtained by using (40) with s set equal 
to r. The special equations developed for this case 
in section 3.2 do not lend themselves well to the 
derivation of the coefficients. 

Let it be emphasized that the p t are the error 
coefficients of the increments of height (Az t ) on the 
layer profile. Of greater interest, perhaps, are the 
error coefficients of the heights (z t ) themselves. 
These coefficients will be called ir h and they are 
defined thus, by analogy with the p z : 



Standard deviation of Zi—o-ki. 



(41) 



The ir i cannot be derived from the p t in any simple 
way because the errors in the calculated Az x are 
mutually correlated even though there is supposed 
to be no correlation between the errors in the virtual 
height data. An expression may be derived for the 
ir h however, similarly to that for the p h using (23) 
as well as (36) and (34); it is 



i}=SJSS v km M nm 

n=r ^U=0m=0 J 



+ ss v km M rm \ y. 

U'=0m=0 J J 



(42) 



Alternatively, in appendix C it is shown that the 
whole analysis can be developed from the start in 
terms of the z t rather than the Az h without altering 
the form of most of the expressions. In the analysis 
so developed, the -k x would be given by the equation 
analogous to (40). 



4. A Test of the Method 

The foregoing analysis has been tested on an 
artificial ionogram. The test consisted of assuming 
a particular layer profile, computing the correspond- 
ing ionogram (fig. 2), and then attempting to recover 
the profile from the data in the ionogram, using first 
the conventional matrix method, and then the vari- 
ous extensions of the method that were developed in 
section 3. In the numerical calculations, the gyro- 
frequency was taken as 1450 kc, and the angle of 
dip as 68.2°, which are the values that apply at 
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Figure 2. The virtual height curves calculated for the assumed 
layer profile. 



Boulder, Colo. The calculations were performed on 
the IBM-650 computer at the Central Radio Pro- 
pagation Laboratory of NBS. 

The assumed layer profile is graphed as the solid 
curve in figures 3 to 5. It is mono tonic, since the 
analysis is restricted to profiles of this type for the 
time being. It represents a distribution of ioniza- 
tion in which there is a ledge of low density beneath 
a layer of high density. Such a distribution might 
arise in the ionosphere soon after dusk, when the 
2?-layer is decaying rapidly but the i^-layer persists. 

The virtual height curves that are derived from 
this profile are shown in the "ionogram" of figure 2. 
The lowest observable frequency (J T ) for the ordinary 
virtual height curve has been taken as 2 Mc, while 
the corresponding lower limit for the extraordinary 
curve is 2.85 Mc. The figure shows only those 
portions of the curves that are supposed to be 
observable. A particularly precise type of numerical 
integration (a 48-point Legendre-Gauss quadrature 
[9]) was used in computing the virtual heights, so 
that these data essentially are free from error. 
Hence any errors that may arise in the deduced layer 
profiles represent systematic defects in the methods 
of analysis. 

The results of various attempts to recover the 
profile are presented in figures 3 to 5. In these 
figures, the points on the deduced profile are shown 
as solid dots, while the broken line is a direct inter- 
polation between adjacent points. 
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Figure 3 shows the results of analyzing the iono- 
gram by the conventional matrix method in which 
only the ordinary virtual height data are used [1, 3]. 
Because the low-density ionization has been ignored, 
all the computed points are much in error, lying 
well above the true profile (the solid curve); the 
error amounts to more than 30 km at the lowest 
observed frequency (2 Mc), while even at 6 Mc it is 
about 6 km. It is this error that the extensions of 
the matrix method seek to eliminate. 
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Figure 1$. The layer profile inferred from analysis of the ordi- 
nary virtual height curve by the conventional matrix method. 

For further explanation, see section 4 of the text. 



However, in attempting to apply the new methods 
various problems were encountered, as follows: 

The first, a comparatively trivial problem, was 
presented by the "rounding-oflf errors" in the digital 
computer. Such errors proved to be particularly 
serious in the inversion of the matrix [Ui m ]. Even- 
tually they were made negligible by performing the 
inversion with a "double-precision" program in 
which the computer carried 18 significant digits of 
each number instead of the usual 8 [10]. 

The main problem proved to be one of balance 
between the effects of random and systematic 
errors on the derived layer profile. The random 
errors are introduced by the inaccuracy of the 
virtual height data, as explained previously. The 
systematic errors are introduced when, in the process 
of converting the original integral equations into the 
corresponding matrix equations, the continuous 
curve of the profile is approximated by a series of 
rectilinear segments (see sec. 2.2.). The need to 
compromise between the two types of error does not 
arise in the analysis of a single virtual height curve, 
that is, in the situation studied in section 2, for in 



that situation the systematic errors can be made 
arbitrarily small by increasing the number of sample 
frequencies and reducing the intervals between 
them. The introduction of a new sample frequency 
adds a new unknown to the set of equations, but 
it also adds a new item of data that is related closely 
to this particular unknown, so that the random 
errors are relatively unaffected. The situation is 
quite different, however, when the virtual height 
curves are incomplete. Then the introduction of a 
new sample frequency in the unobserved range 
increases the number of unknowns without adding 
to the data, so that the precision with which the 
profile can be determined, in the presence of random 
errors, suffers accordingly. It is clear that the 
disposition of the sample frequencies, particularly 
of those in the unobserved frequency range, is a 
delicate problem and one that involves balancing the 
effects of random and systematic errors against each 
other. In view of the complexity of the equations, 
it seemed futile to attempt an analytic solution to 
this problem; instead, the sample frequencies were 
selected by a process of trial and error. The im- 
portance of the random errors was judged from the 
size of the error coefficients, and that of the sys- 
tematic errors from the accuracy with which the 
derived profile reproduced 1 he 1 rue one. An account 
of the series of trials is given now . 

In the search for a satisfactory set of sample 
frequencies, five particular cases were tried in suc- 
cession. Their properties are summarized in table 
1, which gives for each case the values of the para- 
meters r and s, the sample frequency j x that cor- 
responds to each index Z, and the error coefficient 
Pi of each height increment Az t on the layer profiles. 
The profiles themselves were computed only for 
those cases where the error coefficients were reason- 
ably small (i.e., <10), since it was found that the 
other cases developed severe "rounding-off errors" 
in the matrix multiplication (36). 

Cases 1 and 2 

In the first two trials a constant interval of 400 
kc was maintained between sample frequencies. 
Case 1 was exactly-determined, with the number 
of data equal to the number of unknowns, while 
case 2 was overdetermined, with twice as many 
data as unknowns. As would be expected, the error 
coefficients for the overdetermined case were smaller 
than those for the exactly-determined case. In 
both cases, however, the error coefficients for the 
height increments in the observed frequency range 
were excessively large, of the order of several hundred. 
This fact suggests that it is not practical to space 
all the sample frequencies uniformly because random 
errors then have too severe an effect on the lower 
part of the deduced profile. Kather, a relatively 
coarse spacing must be employed in the unobserved 
frequency range, so that in this range the number of 
sample frequencies is reduced and the unknowns 
become very much overdetermined. 
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Table 1. The error coefficients (pi) of the increments of height 
on the layer profile. 

Each coefficient is given in the form of a number, followed in parentheses by 
the power of ten by which that number should be multiplied; these coefficients 
are dimensionless. For further explanation, see section 4 of the text. 



Case 3 

Accordingly, in case 3 the effect was tried of repre- 
senting the entire unobserved frequency range by a 
single interval of sample frequency; in other words, 
the variation of height with plasma frequency was 
assumed to be linear from to 2 Mc. The sample 
frequencies in the unobserved range, however, were 
the same as those in case 2. The error coefficients 
now became reasonably small, and so the layer pro- 
file was computed for this case; the result is shown 
in figure 4. Although this profile is slightly more 
accurate than that of figure 3, it still contains large 
systematic errors which are due to the crudity of 
the straight-line approximation to the layer profile 
in the unobserved frequency range. Evidently a 
representation has to be found for this part of the 
profile that is intermediate in complexity between 
those of cases 2 and 3. 

Case 4 

In case 4, therefore, three sample frequencies 
were placed in the unobserved range, as against 
five for case 2, and one for case 3. Also, the number 
of sample frequencies in the observed range was 
increased, and their spacing was reduced to 200 kc 
in the hope that these changes would help to reduce 
the error coefficients. In fact, however, the error 
coefficients for index and 1 turned out to be dis- 
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Figure 4. The layer profile inferred from both virtual height 
curves by the extended matrix method in case 8. 

For further explanation, see section 4 of the text. 



appointingly large. The conclusion to be drawn from 
this result is that no reliable detail can be obtained 
at the low-frequency end of the unobserved range. 

Case 5 

In case 5, the final trial, the sample frequencies 
were the same as those in case 4, but the ionization 
in the lowest interval of sample frequency was 
ignored: from zero frequency up to j x (0.8 Mc) the 
layer profile was assumed to be horizontal. This 
restriction was imposed on the profile by modifying 

the matrices [M nm ] and [M nw J. The modifications 
consisted of omitting the row and column with 
index zero, and then defining the elements of the 
first column to be equal to unity. These changes 
had the desired effect of reducing the error co- 
efficients to acceptable values, while at the same 
time bringing the inferred layer profile (fig. 5) into 
close agreement with the true profile. The only 
outstanding discrepancy is at the extreme low- 
frequency end of the unobserved range, where the 
assumption that the profile is horizontal makes it 
inevitable. Throughout the observed frequency 
range, the systematic errors of the computed points 
are all less than 1 km. Hence, in this final case, the 
choice of sample frequencies is deemed satisfactory. 

No doubt the method could be improved further 
by continuing this process of trial and error, but the 
results already achieved are felt to give a sufficient 
demonstration of its potentialities. 
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Figure 5. The layer profile inferred by the extended matrix 
method in case 5. 

For further explanation, see section 4 of the text. 



5. Discussion 

This section contains a discussion of one difficulty 
that might be encountered in the practical use of the 
methods of this paper. The difficulty, noted pre- 
viously in section 3.4, is that caused by "lateral 
deviation. " It is well known that, in vertical sound- 
ing, the paths of the ordinary and extraordinary rays 
within the ionosphere are deviated laterally in 
opposite directions [11]. So long as the ionosphere 
is uniform in the horizontal plane, this deviation has 
no effect on the virtual heights. But whenever there 
are horizontal gradients of electron density, then the 
ordinary and extraordinary rays are propagated 
through different distributions of ionization, and for 
neither mode does the virtual height curve correspond 
to the vertical distribution of ionization directly above 
the sounder. This effect has been detected already 
as a source of error in attempts to measure the gyro- 
frequency in the i^-layer by comparison of the pene- 
tration frequencies for the two modes [11]. Pre- 
sumably it also would cause errors in attempts to 
deduce the vertical distribution of ionization by 
joint use of both virtual height curves. 

In this connection, it is of interest that the com- 
parison of the penetration frequencies observed on the 
ionograms indicates whether or not a horizontal 
gradient of electron density is present. If this com- 
parison yields the correct value for the gyrofre- 
quency, then there can be no horizontal gradient; 
such an ionogram could be analyzed with confidence 
to obtain the vertical distribution of ionization. If, 
on the other hand, there is a discrepancy between the 



calculated and true gyro frequencies, then the size of 
this discrepancy provides a measure of the gradient; 
perhaps this information could be used to apply some 
first-order correction to the virtual height curves, so 
as to make them more amenable to analysis. 

6. Conclusions 

In principle, given an ionogram in which the 
virtual heights are not recorded below some lower 
limiting frequency, it is possible to deduce the 
complete layer profile for the full range of frequency 
by combinining the data in the ordinary and extraor- 
dinary virtual height curves; it must be presupposed, 
however, that the profile is monotonic and that the 
data are exact. 

In practice the data contain random errors of 
measurement, and these errors have a severe effect 
on the inferred profile if too much information is 
sought about its shape in the unobserved frequency 
range. To keep the effect of random errors within 
reasonable bounds, detail must be sacrificed in the 
unobserved frequency range, particularly at the 
lowest frequencies; this is to say that the lower part 
of the profile must be characterized by a few param- 
eters only, and that these parameters must be grossly 
overdetermined by the data. The abandonment of 
this low-frequency detail involves systematic errors 
even when the data are exact. In fact, a balance 
has to be struck between random and systematic 
errors. It is possible to strike a balance that gives a 
considerable improvement over the conventional 
matrix method, in which the ordinary virtual-height 
curve alone is analyzed. 

The methods need further extension to cover 
situations where the layer profile is multivalued. 
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7.1. 



7. Appendixes 

Appendix A: Formulas for Computing the Group 
Refractive Index Matrices 



The elements of the matrices [M nm ] and [M nm \ for 
l^m^n, are given by the integrals in (18) and (20), 
respectively. These integrals cannot be evaluated 
analytically, in general, so they have to be evaluated 
numerically. As numerical integrations, however, 
they are difficult to perform for elements on or near 
the diagonal of the matrix (mc^n), since the plasma 
frequency j N then ranges close to the critical value 
for reflection, where the group refractive index 
becomes infinite. Nonetheless, the elements them- 
selves are all finite; and, in fact, each integral may 
be converted into a form that is usable right up to 
the diagonal by a simple change of variable [1]. 
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The procedures for the ordinary and extraordinary 
matrices are slightly different, so they will be treated 
separately. 

a. Ordinary Matrix 

For the ordinary mode, the matrix elements are 
given by the integral 

M nm ^^j- f fm / (UnW* (1 S m S *) . (43) 

To obtain a form in which the integrand stays finite 
as f N approaches J n , the variable of integration is 
changed to 

4>=sm-\f N /f n ). (44) 

Then 

M nm =^~ f* m '{(l-3T»)V(/»,/»sin $}d<t> 

(lgmgw), (45) 

where <$> m , 3> m _i are the values of <t> for which f N =f m , 
f m -i, and where 



X n =(f N /f n y=sm% 



(46) 



If the frequency interval Af m is a constant, say A/, 
so ihditj n =n A/, then 



M« 



%f*' {(1-XJV (/»,/» sin *)}<Z0 

(l^m^n). (47) 



In evaluating this integral, the integrand is com- 
puted from the equations of the magneto-ionic 
theory, using the familiar variables 



Y=M 



Y L =Y cosd Y T =Ysmd 



(48) 
(49) 



in which 6 is the angle between the direction of 
propagation (the vertical) and the earth's magnetic 
field; the effect of collisions is ignored, as is usual in 
this type of work. From these variables, the 
following parameters are defined : 5 



S={Y* T +iYl(l-Xy}i, 
D=2(l-X)-E, 
E=rY£-kS, 
F=E+2Yi{l-X). 



(50) 
(51) 
(52) 
(53) 



* These parameters are similar to a set introduced by Shinn [12] who is respon- 
sible for several of the formulas quoted here. However, some of Shinn's 
definitions have been altered slightly. 



It is convenient to compute the integrand as the 
product of two factors : 



(l-Z)V={w 



'>m s 



(54) 



Here /x is the wave refractive index for the ordinary 
mode. Each of these factors, separately, also 
remains finite as/^ approaches /. They are given by 



MM 



and 



'=M 2 +J{ 



E(E+2YD- 



4Fi(l-Z)(l-Z 2 ) 



S 



l-X 



= 1- 



while, furthermore 



,2=1 



2YJX 
'E+2YI' 



EX 



} 

(55) 
(56) 

(57) 



b. Extraordinary Matrix 



The integral for the elements of the extraordinary 
matrix is 



M nm =- 



=-7g- m f^xifn, 
Af m Jf m-l 



Jn) d/N 



(l^m^n) (58) 



where 



fn = ^fH[{l+MfJfH)} 2 ±l] 



(59) 



and the sign is to be taken positive or negative 
according to whether f n is greater or less than j H . 
Equation (58) may be changed into a form suitable 
for numerical integration by making the substitu- 
tion 



Then 



J J *— 1 



1 

A/, 



4>=sm- 1 (f N /f n ). 



(1— X n ) *m« (/»,/» sin 4>)}d<t> 



(60) 



(lgmgri) (61) 



where <t> m ,4> m -i are the values of <j> for which fit=l m , 
j m -\, and where 



X n =(fN/Ly=sm 2 4>. 



(62) 



For the case where the frequency interval Af m is a 
constant, say A/, 



M n 



J*m-l 



(1-X n ^u.' x (f n ,j n sm<j>)}d^ (lSf»£n). 

(63) 
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Again, the integrand is computed as the product 
of two factors. Here the factors are, firstly, 



Mx»x = vl+jj\ i— A— X 2 

and secondly, the square root of 

l-X D 

2(1TY)(1-X±T) 



(l-x 2 ) (i-/4) it 



s 



] 



(64) 



A 



(1 + 2YI/E). (65) 



The wave refractive index for the extraordinary 
mode is given by 



/4 = 



2XQ.-X) 
D 



Note also that 



X=X(ITY). 



(66) 



(67) 



Here, and in (65) above, the upper sign applies to 
frequencies above f H , and the lower to frequencies 
~belowf H . 

7.2. Appendix B: Treatment of the Overdetermined 
Case by the Method of Least Squares 



The equations to be solved are 



K^J2 M nl Az l (r<n<r+s-l) 

1=0 



for the ordinary mode, and 



h n c^P 1 M nl Az l 

1 = 



(r<n<r-{-s—l) 



(17) 



(19) 



for the extraordinary. Observed values of the 

h' n and h n are given for r^n^r-\-s— 1, where s^r, 
but these data contain random errors of measure- 
ment The aim is to find a set of Az t which, if 
inserted in the equations, would yield calculated 

values for the h' n and h n that agree with the observed 
values as closely as possible. The measure of dis- 
agreement is the sum of the squares of the differences 
between the two sets of values, and this quantity 
is to be minimized. 

Now, if certain values are assumed for the Azi t 
the implied errors of measurement are 



and 



1 = 



E n =l n -it l M nl Az l . 

1=0 



(68) 



(69) 



For the ordinary mode, the sum of the squares of the 
implied errors is 

r+s— 1 r -\-s— 1 r-j-s— 1 n 

S E*= s (KV-2 2 IbMmhfa, 

n=r n = r n=r 1=0 

+ r £ 1 SEM.AMA,. (70) 

n = r k = 1=0 

A similar expression applies to the extraordinary 
mode. These two expressions must be added, and 
the total minimized. The condition thai it shall be 
a minimum with respect to variations of all the 
Az m is that 



d(A2 : 



r+s-l „ 



(71) 



m) n = r 



for all m (0 ^m^r+6'— 1). On evaluating the 
derivatives, this condition yields the Az t as the 
solutions of (33). 

7.3. Appendix C: The Analysis Presented in Terms 
of the Ordinates z, 

In the analysis described in this paper, the ordi- 
nates of the layer profile have been arrived at by 
two distinct steps: First, various matrix equations 
have been set up and solved in terms of the height 
increments Az h then these increments have been 
summed to obtain the actual ordinates z x . It is 
possible to combine these two steps by rewriting the 
matrix equations with the ordinates themselves as 
the unknowns. Consider, for instance, the eq (17) 
for the ordinary virtual height. It may be re- 
arranged as follows: 



K^ZM nm Az m 



±M n0 AZ +JlM nm AZn 



~M n0 z + S M nm (z m — z m - 2 ) 

m=l 



iM nO z +^M nm z m —^l M nm z m -i 

m = \ m = 1 



(17) 





— S i -M-nrnZm. / i M-n, m.4--\Z m 

m=Q m=0 




n 
m = 


say, 




where 






(M nm —M n<m+1 (m<n) 




o n m=^M nn (m=n) 




to (m>n). 



(72) 



(73) 
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This is the form in which Budden [1] uses the matrix 
equation; it is, of course, completely equivalent to 
the other. 

The eq (19), for the extraordinary virtual height, 
may be rearranged in the same way: 



l^n — / i ^nm^n 



say. 



(74) 



Since the eqs (72) and (74) are of the same form as 
(17) and (19), it is clear that the analysis of this paper 
could have been based equally well on the new equa- 
tions, and developed in terms of the ordinates z%\ the 
arguments would proceed just as before, with the z x 
replacing the Az h so there is no need to give the details. 
The choice of which of the two sets of equations to 
start from is largely a matter of convenience, though 
it may be that, in computation, one would involve 
smaller "rounding-off errors" than the other. 
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